Simulation of merging binary neutron stars in full general relativity: F = 2 case 
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We have performed 3D numerical simulations for merger of equal mass binary neutron stars in full 
general relativity. We adopt a F-law equation of state in the form P = (F — l)pe where P, p, 
e and T are the pressure, rest mass density, specific internal energy, and the adiabatic constant 
with r = 2. As initial conditions, we adopt models of corotational and irrotational binary neutron 
stars in a quasi-equilibrium state which are obtained using the conformal fiatness approximation 
for the three geometry as well as an assumption that a helicoidal Killing vector exists. In this 
paper, we pay particular attention to the final product of the coalescence. We find that the final 
product depends sensitively on the initial compactness parameter of the neutron stars : In a merger 
between sufficiently compact neutron stars, a black hole is formed in a dynamical timescale. As 
the compactness is decreased, the formation timescale becomes longer and longer. It is also found 
that a differentially rotating massive neutron star is formed instead of a black hole for less compact 
binary cases, in which the rest mass of each star is less than 70 — 80% of the maximum allowed 
mass of a spherical star. In the case of black hole formation, we roughly evaluate the mass of the 
disk around the black hole. For the merger of corotational binaries, a disk of mass ~ 0.05 — O.IM, 
may be formed, where M* is the total rest mass of the system. On the other hand, for the merger 
of irrotational binaries, the disk mass appears to be very small : < 0.01M« . 
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I. INTRODUCTION 

Several neutron star-neutron star binaries are known 
to exist in our galaxy According to high precision 
measurements of the general relativistic (GR) effects on 
their orbital motion, three of these binaries are going 
to merge as a result of gravitational radiation reaction 
within the Hubble timescale ~ 10^° years. What are 
the final fates of these binaries after the mergers? Their 
total gravitational masses are is in a narrow range ~ 
2.65 — 2.85M0 where Mq denotes the solar mass. The 
stars will not be tidally disrupted before the merger, since 
the masses of the two stars in each binary are nearly 
equal. Hence, the mass loss from the binary systems is 
expected to be small during their evolution and the mass 
of the merged object will be approximately equal to the 
initial mass. The maximum allowed gravitational mass 
for spherical neutron stars is in a range ~ 1.5 — 2.3M0 
depending on the nuclear equation of state ||^,^. Even 



if we take into account the effect of rigid rotation, it is 
increased by at most ~ 20% [Q. Judging from these facts, 
the compact objects formed just after the merger of these 
binary systems seem bound to collapse to a black hole. 

However, this is not the case if the merged object ro- 
tates differentially. The maximum allowed mass can be 
increased by a larger factor (> 50%) due to the differen- 
tial rotation which suggests that the merged objects 
of ~ 2.5 — 3M0 may be dynamically stable against grav- 
itational collapse to a black hole. Such differentially ro- 
tating stars could be secularly unstable, since viscosity or 
magnetic field could change the differential rotation into 
rigid rotation. A star with a high ratio of rotational en- 
ergy to the gravitational binding energy could also be sec- 
ularly unstable to gravitational wave emission Q . These 
processes might dissipate or redistribute the angular mo- 
mentum, and induce eventual gravitational collapse to a 
black hole. However, the timescales for such secular in- 
stabilities are in general much longer than the dynamical 
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timescale of the system. Hence, the merged objects may 
not collapse to a black hole promptly, but remain as a 
massive neutron star supported by differential rotation 
at least for these secular timescales. These facts imply 
that the final product of the merger of binary neutron 
stars is an open question depending not only on the nu- 
clear equation of state for high density neutron matter 
but also on the rotational profile of the merged object. 

Interest in the final product of binary coalescence has 
been stimulated by the prospect of future observation of 
extragalactic binary neutron stars by gravitational wave 
detectors. A statistical study shows that mergers of bi- 
nary neutron stars may occur at a few events per year 
within a distance of a few hundred Mpc . This suggests 
that binary merger is a promising source of gravitational 
waves. Although the frequency of gravitational waves in 
the merging regime will be larger than IkHz and lies be- 
yond the upper end of the frequency range accessible to 
laser interferometers such as LIGO |^ for a typical event 
at a distance ~ 200Mpc, it may be observed using spe- 
cially designed narrow band interferometers or resonant- 
mass detectors |9| . Such future observations will provide 
valuable information about the merger mechanism of bi- 
nary neutron stars and the final products. 

Interest has also been stimulated by a hypothesis about 
the central engine of 7-ray bursts (GRBs) [0. Recently, 
some of GRBs have been shown to be of cosmological ori- 
gin In cosmological GRBs, the central sources must 
supply a large amount of the energy > 10^^ ergs in a very 
short timescale of order from a millisecond to a second. 
It has been suggested that the merger of binary neutron 
stars is a likely candidate for the powerful central source 
[0. In the typical hypothetical scenario, the final prod- 
uct should be a system composed of a rotating black hole 
surrounded by a massive disk of mass > O.IMq, which 
could supply the large amount of energy by neutrino pro- 
cesses or by extracting the rotational energy of the black 
hole. 

To investigate the final product of the merger theo- 
retically, numerical simulation appears to be the unique 
promising approach. Considerable effort has been made 



for this in the framework of Newtonian and post- 
Newtonian gravity ||l^-[l7|]. Although these simulations 
have clarified a wide variety of physical features which are 
important during the coalescence of binary neutron stars, 
a fully GR treatment is obviously necessary for determin- 
ing the final product because GR effects are crucial. 

Intense efforts have been made for constructing a reli- 
able code for 3D hydrodynamic simulation in full general 
relativity in the past decade [p^^. Recently, Shibata 
presented a wide variety of numerical results of test prob- 
lems for his fully GR code and showed that simulations 
for many interesting problems are now feasible pl| . 

To perform a realistic simulation, we also need re- 
alistic initial conditions for the merger, i.e., a realistic 
density distribution and velocity field for the component 
stars. Since the timescale of gravitational wave emission 
is longer than the orbital period, we may consider the 
binary neutron stars to be in a quasi-equilibrium state 
even just before the merger. The velocity field in the 
neutron stars has been turned out to be nearly irrota- 
tional because (1) the shear viscosity is too small to re- 
distribute angular momentum to produce a corotational 
velocity field during the timescale of gravitational wave 
emission and (2) the initial vorticity of each star is neg- 
ligible as long as the rotational period of the neutron 
stars is not ~ milliseconds [Q. Therefore, as realistic 
initial conditions, we should prepare a quasi-equilibrium 
state of binary neutron stars with an irrotational veloc- 
ity field. Recently, several groups have developed various 
numerical methods to compute GR irrotational binary 
neutron stars in quasi-equilibrium in a framework with 
the appropriate approximations that a helicoidal Killing 
vector (see Eq. ( ^.l| ) for definition) exists and that the 
three geometry is conformally flat [p3|-p5| . Their numer- 
ical results have been compared and found to agree to 
within a few percent error in the gravitational mass and 
the central density as a function of orbital separation. 

In this paper, we perform simulations for the merger 
of binary neutron stars of equal mass by using these new 
numerical implementations developed recently. As a first 
step, we adopt a simple F-law equation of state with 
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r = 2 as a reasonable qualitative approximation to the 
high density nuclear equation of state. Although micro- 
scopic effects such as cooling due to neutrino emission or 
heating due to bulk viscosity may become important for 
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discussing the merging process in detail |16|, we neglect 
them here for simplicity. The purpose of this paper is to 
investigate the dynamical nature of the mergers, the final 
products after the mergers and the dependence of these 
outcomes on the initial velocity field and the compact- 
ness parameter of the binary neutron stars. Simulations 
are performed not only for irrotational binaries but also 
for corotational ones to clarify the difference due to the 
initial velocity field. 

Throughout this paper, we adopt the units G = c = 
Mq — 1 where G and c denote the gravitational constant 
and speed of light, respectively. Latin and Greek indices 
denote spatial components (1 — 3) and space-time com- 
ponents (0 — 3), respectively. 6^^) denotes the Kro- 
necker delta. We use Cartesian coordinates x'^ = {x, y, z) 
as the spatial coordinates with r ~ -^/x^ + + z^; t de- 
notes coordinate time. 

II. METHODS 
A. Summary of formulation 

We have performed numerical simulations using the 
same formulation as in pH , to which the reader may re- 
fer for details about the basic equations, the gauge condi- 
tions and the computational method. The fundamental 
variables used in this paper are: 

p: rest mass density, 
e: specific internal energy, 
P: pressure, 
u^: four velocity, 

" = — ' 

a: lapse function, 
p^: shift vector, 

7y : metric in 3D spatial hypersurface. 

Ho) 



Kij: extrinsic curvature. 

These variables (together with auxiliary functions Ei and 
the trace of the extrinsic curvature KfJ^) are numerically 



evolved as an initial value problem (see |26 2^] for de- 
tails of the numerical method for handling the evolution 
equations and initial value equations). Several test cal- 
culations, including spherical collapse of dust, stability 
of spherical neutron stars, and the evolution of rotating 
neutron stars as well as corotational binary systems have 
been presented in Violations of the Hamiltonian 

constraint | |28[ | and conservation of rest mass and angu- 
lar momentum pof are monitored to check the accuracy, 
and we stop simulations before the accuracy becomes too 
poor. Black holes which are formed during the last phase 
of the merger are located with an apparent horizon finder 
described in |30| . 

We also define a density p*(= pau^e^'^) from which the 
total rest mass of the system can be integrated as 



= / (f'xp: 



(2.1) 

We have performed the simulations using a fixed uniform 
grid with a typical size of 233 x 233 x 117 for the x — y — z 
directions, respectively, and assuming refiection symme- 
try with respect to the z — Q plane. All of the results 
shown in Sec. IV are obtained from simulations with this 
grid size. We have also performed a number of test simu- 
lations with a smaller grid size of 193 x 193 x 97 changing 
the grid spacing and location of the outer boundaries to 
confirm that the results do not change significantly. For 
one model (12) (see Sec. Ill and Table I), we performed 
a simulation with a larger grid size of 293 x 293 x 147, 
widening the computational domain to investigate the 
effect of the outer boundary on the gravitational wave- 
forms. 

The slicing and spatial gauge conditions which we use 
in this paper are basically the same as those adopted 



7= det(7y) = e^^-A = V'', 



in 1 26 27|pl[ |; i.e., we impose an "approximate" maximal 
slicing condition {K^ 0) and an "approximate" mini- 
mum distortion gauge condition {Di{dtj^^) — where Di 
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is the covariant derivative with respect to 7ij). However, 
for the cases when a merged object cohapses to form a 
black hole, we slightly modify the spatial gauge condi- 
tion in order to improve the spatial resolution around 
the black hole forming region. The method of the modi- 
fication is described in Sec. II.B. 

Throughout this paper, we assume a F-law equation of 
state in the form 



P = (r - l)pe, 



(2.2) 



where F is the adiabatic constant. For the hydrostatic 
problem, which appears in solving for initial value con- 
figurations, the equation of state can be rewritten in the 
polytropic form 

P = Kp^, F = l + -, (2.3) 
n 

where X is a constant (different from Kf,'') and n is the 
polytropic index. We adopt F = 2 (n 1) as a qualita- 
tive approximation to realistic, moderately stiff equations 
of state for neutron stars. 

Physical units enter only through the constant K, 
which can be chosen arbitrarily or completely scaled out 
of the problem. In the following, we quote values for K = 
200/7r, for which in our units (G = c = Mq = 1) the ra- 
dius of a spherical star is i? = {nK/2y/^ — 10(^ 15km) 
in the Newtonian limit. Since if"/^ has units of length, 
dimensionless variables can be constructed as 

J = JK--,Porh = PorhK-^'^P = pK", (2.4) 

where Mg, J and Poib denote the gravitational mass, an- 
gular momentum and orbital period. All results can be 



rescaled for arbitrary K using Eqs. (2.4). (For example, 
the maximum value of Af* for spherical stars is 1.435M0 
for K = 200/n{G^Ml/c^) with the rest mass density 
~ 3.1 X lO^^g/cm''. If we want to choose Af* as 2Mq, 
we should change K to 123.7(G'^Af|,/c^) and the corre- 
sponding density is then ~ 1.6 x lO^^g/cm'^; cf. Fig. 1. ) 
In other words, the invariant quantities are only the di- 
mensionless quantities such as A'h/Mg, Mg/R, Mg/Pah 
and J/A//2. 



B. Spatial gauge condition 

When a black hole is not formed as a final product 
of the merger, we adopt the approximate minimum dis- 
tortion gauge condition as our spatial gauge condition 
(henceforth referred to as the AMD gauge condition). 
However, as pointed out in previous papers ||2l| , p7|l , dur- 
ing black hole formation (i.e., for an infalling radial ve- 
locity field), the expansion of the shift vector 9;/?" and 
dtcj) — di(3^/6 become positive using this gauge condition 
together with our slicing Kf.'' ~ 0. Accordingly, the coor- 
dinates diverge outwards and the resolution around the 
black hole forming region becomes worse and worse dur- 
ing the collapse. This undesirable property motivates us 
to modify the AMD gauge condition when we treat black 
hole formation. Following |^ , we modify the AMD shift 
vector as 



P' = PXMB-fit,r)—-PlMD- (2.5) 
r + e 

Here /?amd denotes the shift vector for the AMD gauge 
condition, /3amd = (^amd / (''" + e is a small con- 
stant much smaller than the grid spacing, and /(i, r) is 
a function chosen specifically as 

1 



/(t,r) = /o(t)- 



(2.6) 



' 1 + {r/3Mgof ■ 

where Mgo denotes the gravitational mass of a system at 

t = Q. We determine fo{t) from (j>o — (j>{r ~ 0). Taking 

into account the fact that the resolution around r = 

deteriorates when (jjo becomes large, we choose /o as 

f 1 for 00 > 0.8, 

fo{t) = <^ 2.500 - 1 for 0.4 < 0o < 0.8, (2.7) 
[ for 00 < 0.4, 

Note that for spherical collapse with /o = 1, 9^/3' ~ 
and 9(0 ~ at r = 0. We employ this modified gauge 
condition whenever a merged object collapses to form a 
black hole. 

It is worth mentioning that with this modification, the 
coordinate radius of the apparent horizon (when it is 
formed) becomes larger than without the modification. 
This implies that more grid points are located along the 
radius of the apparent horizon and accuracy for determi- 
nation of the apparent horizon is improved. 
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III. INITIAL CONDITIONS 

Even just before the merger, the binary neutron stars 
are considered to be in a quasi-equihbrium state be- 
cause the timescale of gravitational radiation reaction 
- 5/{64f^(A/grj)5/3| where n denotes the orbital an- 
gular velocity of the binary neutron stars, is several times 
longer than the orbital period. Thus, for performing a 
realistic simulation of the merger, we should prepare a 
quasi-equilibrium state as the initial condition. In this 
paper, we construct such initial conditions as follows. 

First, we assume the existence of a helicoidal Killing 
vector 



Lii — 



(3.1) 



Since emission of gravitational waves violates the heli- 
coidal symmetry, this assumption does not strictly hold 
in reality. However, as mentioned above, the emission 
timescale of gravitational waves is several times longer 
than the orbital period even just before the merger (c/. 
Table I) so that this assumption can be acceptable for 
obtaining an approximate quasi-equilibrium state. In ad- 
dition to this assumption, we adopt the so-called confor- 
mal flatness approximation in which the three geometry 
is assumed to be conformally flat, for simplicity. 

In this paper, we consider irrotational and corota- 
tional binary neutron stars. Then, the geometric | |3^ ] 
and hydrostatic equations p4| for solutions of the quasi- 
equilibrium states are described as 



= -27r(pW - Kp^)i,^ - i^S'^'S^'U^Lku (3.2) 

o 



5,, A/3' + - 2L^,S'''(d,a - ^d,^ 



ah 
w 

where 



3'- —J- V 
— IGnaphuJUj, 

hukV^ = const.. 



h = i + Krp^-^/{r-i), 



(3.3) 

(3.4) 
(3.5) 



(3.6) 
(3.7) 



2a 



yk ^_pk 



kl ""i 



t. 



(3.8) 
(3.9) 



and A denotes the flat Laplacian in the three space. Eqs. 
(^^-^^4) are the geometric equations and Eq. (3_^) is the 
so-called Bernoulli equation. can be regarded as the 
coordinate three velocity in the corotating frame rotating 
with angular velocity f2. 

In the case of corotational binaries in which V*' = 0, 
Ui is written as 

u, = w^p\e,,k^x'' + S,,l3^)/a. (3.10) 

In the case of irrotational binaries, on the other hand, Ui 
is written as 



(3.11) 



where $ denotes the velocity potential which satisfles an 
eUiptic PDE @ 

S'^d,{pai;'^h-^dj^) - d,[pah^^ [t + /?')] = 0, (3.12) 

with the following boundary condition at the stellar sur- 
face; 



= 0. 



(3.13) 



The above Poisson type equations such as Eqs. ( 3.2 



3.5) and (3.12) as well as the Bernoulli equation (3.5) are 
solved iteratively with appropriate boundary conditions. 
Corotational binaries are calculated using the same nu- 
merical method as adopted in pl|] . Irrotational binaries 
are calculated using the method developed recently by 
Uryii and Eriguchi ( ||2^ ], to which the reader may refer 
for details). 

For the corotational case, we prepare binaries with sev- 
eral compactness parameters, with the surfaces of the two 
stars coming into contact. As shown in such bina- 
ries with r = 2 are located near to the energy minimum 
along the sequence of corotational binaries of constant 
rest mass. Therefore, they are expected to be located 
near to a marginally stable point for hydrodynamic p4| ] 
or GR orbital instability. 
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For the irrotational case, the sequence of binaries of 
constant rest mass ends when cusps (i.e., Lagrange (LI) 
points) appear at the inner edge of the stars ||2^,Q . This 
is the case for any compactness parameter. If the stars 
in the binary system approach further, mass transfer will 
begin and the resulting state is not clear. As shown in 
, the closest binaries with cusps are far outside the 
energy minimum for F = 2, which indicates that they 
are stable against hydrodynamic and GR orbital insta- 
bility. Furthermore, the gravitational radiation reaction 
timescale is several times longer than the orbital period 
(c/. Table I). Thus, if we choose such a binary as the 
initial condition for a simulation, a few orbits are main- 
tained stably before the merger starts, decreasing the 
orbital separation and changing the shapes in a quasi- 
adiabatic manner. 

It is still difficult to perform an accurate simulation for 
such a quasi-adiabatic phase. It is desirable to choose 
a binary state which is located near to the unstable 
point against hydrodynamic or GR orbital instability and 
starts merging soon; i.e., a state after the nearly adia- 
batic phase. However, a method for obtaining such a 
state has not yet been developed. Hence, in this paper, 
we prepare the following initial conditions modifying the 
quasi-equilibrium state slightly. First, we prepare a bi- 
nary in which cusps appear at the surfaces. Then, we 
reduce the angular momentum by ~ 2.5% from the quasi- 
equilibrium state to destabilize the orbit and to induce 
the merger promptly . We deduce that such an initial 
condition can be acceptable for the investigation of the 
final products after the merger since the decrease factor 
is still small. We performed test simulations changing 
the decrease factor slightly in the range 2 — 3%, and we 
indeed found that the results shown in Sec. IV are only 
weakly dependent on this parameter. 

In the numerical computation of the quasi-equilibrium 
states as initial conditions, we typically adopt a grid spac- 
ing in which the major diameter of each star is covered by 
~ 35 grid points. With this resolution, the error for the 
estimation of Alg and J is less than 1% (e.g., |Q), and 
the gravitational radius of the system defined as GMgo / 



is covered with ^4 — 7 grid points (c/. Table I). Keeping 
this grid spacing, the outer boundaries of the computa- 
tional domain in the simulation are located at ^ O.SAgw 
with 233 X 233 x 117 grid points, where Agw(= tt/H) de- 
notes the characteristic wavelength of gravitational waves 
emitted from the binaries in a quasi-equilibrium state. 
With this setting, gravitational waveforms are not ac- 
curately evaluated because the outer boundaries are not 
located in the wave zone. In this paper, we pay par- 
ticular attention to the merger process, final products, 
and dependence of these outcomes on initial parameters 
of the binaries, but we do not treat the accurate extrac- 
tion of gravitational waveforms and hence the accurate 
computation of gravitational radiation back reaction. As 
mentioned above, we start with binaries in almost dy- 
namically unstable orbits. This implies that the effect 
of radiation reaction is not very important in the early 
phase of the merger. In the later phase of the merger, the 
dynamical timescale seems to be shorter than the emis- 
sion timescale of gravitational waves and the evolution of 
the merged object due to gravitational radiation is secu- 
lar. Hence, we deduce that the effect due to the error in 
evaluating the radiation reaction is small throughout the 
evolution. 

In Fig. 1, we show the relation between the rest mass 
and the maximum density pmax of each star for bi- 
naries in quasi-equilibrium states. The solid line denotes 
the relations for spherical neutron stars. The crosses and 
filled circles denote those for the corotational and irro- 
tational binary neutron stars, respectively. The binaries 
which are used in the following simulation as initial con- 
ditions are marked with (CI), (C2), (C3), (II), (12) and 
(13) (C and I denote "corotational" and "irrotational", 
respectively). The relevant quantities for these initial 
conditions are shown in Table I. We note that the orbital 
period is calculated as 

Porb = 1.5msec ( ( . (3.14) 

In a realistic situation, each star of the binary has a 
small approaching velocity because of gravitational radi- 
ation reaction. We approximately add this to the above 
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quasi-equilibrium state in setting the initial conditions. 
According to the quadrupole formula with Newtonian 
equations of motion, the absolute value of the approach- 
ing velocity of each star is written as Q 

Wa^ 1. 6(1/^0^^)2 . (3.15) 

Thus, in giving initial conditions, we change Ui to be 

\x\ 

Ux ^ {Ux)cq_- Va. — , (3.16) 

where {ux)cci denotes of the quasi-equilibrium state. 
Here, we implicitly assume that the center of mass of 
each star is initially located on the x-axis (see Figs. 2-4 
and 9-11). 

For models (CI), (C2) and (C3), we performed simu- 
lations without the approaching velocity and found that 
the outcomes such as the final products and the disk mass 
depend only weakly on this approaching velocity. Thus, 
it does not seem to affect the following results signifi- 
cantly. 

IV. NUMERICAL RESULTS 
A. corotational cases 

In Figs. 2-4, we show snapshots of the density contour 
lines for and velocity field {v^ ^v"^) in the equatorial 
plane at selected times for models (CI), (C2), and (C3), 
respectively. For (CI), a new massive neutron star is 
formed, while for other black hole is formed. We 

note that for model (C2), we could not determine the 
location of the apparent horizon before the simulation 
crashed, because the grid spacing was too wide to satis- 
factorily resolve the black hole forming region. However, 
the central value of the lapse function is small enough 
< 0.01 at the crash so that we may judge that a black 
hole is formed in this simulation. On the other hand, for 
model (C3), we can determine the location of the appar- 
ent horizon (see the thick solid line in the last snapshot 
of Fig. 4). 

Irrespective of the compactness parameters, the orbital 
distance gradually decreases due to the initial approach 



velocity. When it becomes small enough to destabilize 
the orbit due to the hydrodynamic or GR orbital in- 
stability, the orbital distance begins to quickly decrease 
and in the outer part, spiral arms are formed. For more 
compact binaries, the decrease rate of the orbital sepa- 
ration is larger because the initial approach velocity is 
larger, and the orbit soon becomes unstable. We deduce 
that the neutron stars for model (C3) are initially lo- 
cated near to the innermost stable circular orbit against 
GR orbital instability because their initial compactness 
Ci = {Mg^n)'^/'^ ~ Mgo/a, where a denotes the orbital 
separation, is nearly equal to 1/6 (see Table I). Indeed, 
they begin merging soon after the simulation is started. 
Once merger begins, the spiral arms continue to develop 
transporting angular momentum outward in the outer 
part of the merged object. 

For model (CI), the inner part first contracts after the 
orbit becomes unstable, but subsequently it bounces due 
to the pressure and centrifugal force. The shape of the 
merged object changes from ellipsoidal to spheroidal, re- 
distributing the angular momentum as well as dissipating 
it by gravitational radiation. Eventually, it forms a new 
rapidly rotating neutron star. In Figs. 5 and 6, we show 
the density contour lines for in the x — z plane and 
the angular velocity = {xv^ — yv^) / (a;^ -I- y^) along the 
X and y-axes in the equatorial plane at t — 2.07Porb- 
It is found that the new neutron star is highly flat- 
tened and differentially rotating HJ]. Note that the 
mass inside r ~ 7.5Mgo which appears to constitute the 
merged object is ~ 0.97Af,(~ 2.16). Since the maxi- 
mum allowed mass of a spherical star with K — 200/7r 
is M^^niax — 1.435, the mass of the new neutron star 
is ^ 50% larger than this |3^]. We point out that we 
have monitored the evolution of K'{x^^) = P/ which 
is initially equal to K anywhere in the star and can be 
regarded as a measure of the entropy distribution. Since 
shock heating is not very effective in the merging, we have 
found that the value of K' increases by at most ~ 10% in 
the regions of high density. (Note that in the low density 
regions such as near to the surface of the merged object, 
K' is slightly larger.) Thus, the role of thermal energy 



7 



increase is not significant for supporting the large mass 
in contrast with the case of head-on collision |36y2C| ] . The 
effect of differential rotation is important in the present 
case. 

We note that the new rotating neutron star has a non- 
axisymmetric structure at the time when we stopped the 
simulation. Therefore it will evolve further as a result of 
gravitational wave emission, and may become unstable 
against gravitational collapse to become a black hole after 
a substantial amount of angular momentum is carried 
away 

For models (C2) and (C3), after the orbit becomes un- 
stable, the inner part contracts due to self-gravity with- 
out bouncing because the pressure and centrifugal force 
are not strong enough to balance the self-gravity. Sub- 
sequently it collapses to form a black hole. Since the 
compactness is sufficiently large for model (C3), the in- 
ner part quickly collapses to form a black hole without 
a significant bounce. On the other hand, in the case of 
model (C2), the formation timescale of the black hole is 
longer because the compactness is smaller. 

To show the features of the collapse around the central 
region, we show a at r = as a function of t/Porb in Fig. 
7. For model (C3), a{r = 0) quickly approaches to zero, 
but for model (C2), the decrease rate becomes small at 
t ~ 1.2Porb- This difference indicates that the collapse is 
decelerated by the pressure and/or centrifugal force. 

We note that for model (C2), the initial value of J/M^ 
is larger than unity. Nevertheless, a black hole appears 
to be formed after the merger. This indicates that some 
mechanisms for angular momentum transfer or dissipa- 
tion act to decrease J/M^ to less than unity during the 
merger. We can expect that the following two mecha- 
nisms are effective, (a) In the case of corotational bina- 
ries, the outer part has a large amount of the angular mo- 
mentum and spreads outwards forming the spiral arms. 
As a result, the specific angular momentum in the inner 
part which finally forms a black hole is smaller than that 
of the outer part and J/Mg can be smaller than unity in 
the inner region, (b) The effect of gravitational radiation 
can reduce the magnitude of J/M^ which is estimated as 



follows: If the system has a characteristic angular veloc- 
ity fie, the relation between the energy loss 6E{> 0) and 
the angular momentum loss SJ{> 0) due to gravitational 
radiation can be written as Qc6J — SE. If we assume 
SE ^ Algo and SJ <C Jo where Jq denotes the initial 
value of J, the resulting J/Mg becomes 



Jo ^ SJ 
{Mgo - SE)^ 
5J_ 

' Jo 



Jo 



1 - 



dJ_ 

Jn 



25E\ 
Wo) 



-1 



/2Jo 



\M'f 



(4.1) 



Here, Jo/MgQ ^ 1, and because of the fact that gravita- 
tional waves are efficiently emitted in the early phase of 
merger, we may set flc ^ ^ and consequently QcMgo <C 1 
(see Table I). Thus, 2Jailc/Mgo (the second term in { } 

ih is ap- 



of Eq. (4.1)) is much less than unity, and Eq 
proximately (Jo/MgQ)(l — SJ/Jq). Using the quadrupole 
formula and the Newtonian expression for the angular 
momentum, SJ/Jq in one orbital period for a binary sys- 
tem of point masses is M 



6J_ 
Jo 



^(AJ,of^)V3 = o.0876(^ 



5/2 



(4.2) 



Therefore, J/Mg can decrease by '-^ 10%. 

Since the gradient of the metric becomes very steep in 
the high density region of the merged object, the simu- 
lation could not be accurately continued for models (C2) 
and (C3) after a at r = becomes less than ^ 10^^. Al- 
though we cannot strictly calculate the final states of the 
disks around the black holes for these models, we may 
extrapolate the final state from the evolution of the cen- 
tral region as follows. In Fig. 8, we show time evolution 
of the fraction of the rest mass inside a coordinate radius 
r, defined as 

M,{r) _ _J_ 



(4.3) 



'\<r 



for models (C2) and (C3). We choose r = 1.5,3 and 
A.5MgQ as coordinate radii. It is found that more than 
95% of the total rest mass is inside r ~ A.bMgQ, and 
a small fraction < 3 — 5% of the total rest mass can 
be in a disk around the black hole at r > 4.5Afgo. For 
model (C3), the location of the apparent horizon is at r ^ 
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1.2Mgo at t = l.OSPorb, so that most of the matter inside 
r = 1.5Mgo seems to be swallowed by the black hole 
eventually. On the other hand, since the newly formed 
black hole seems to be rotating rapidly {J/Mg ~ 0.8—0.9, 
see Table I) , the innermost stable circular orbit is located 
near the event horizon and so even some of the matter 
located between r = l.5MgQ and 3MgQ may go to form 
the disk around the black hole. Hence, it may still be 
possible that a very compact disk of mass ^ 0.05M* and 
radius ~ 3Afgo is formed eventually for model (C3). For 
model (C2), we cannot discuss details because we could 
not determine the apparent horizon. As shown in Fig. 8, 
the fraction of matter inside r — l.dMgo is still increasing 
at the time when we terminated the simulation, so that 
the mass fraction of the compact disk at r 3Mgo seems 
to be at most 0.05Af*. Thus, a disk of mass at most 
~ 0.05 — O.IA/* may be formed around black holes in an 
optimistic estimation. 

B. irrotational cases 

In Figs. 9-11, we show snapshots of the density con- 
tour lines for and velocity field (v^jV^) in the equa- 
torial plane at selected times for models (II), (12), and 
(13), respectively. For model (II), a new massive neutron 
star is formed, while for the other black hole is 

formed. We note that we could not determine the loca- 
tion of the apparent horizon for model (12) before the 
simulation crashed. However, the central value of the 
lapse function is small enough < 0.01 at the crash, so 
that we judge that a black hole is formed in this simu- 
lation as in the case (C2). On the other hand, we could 
determine the location of the apparent horizon for model 
(13). 

As in the corotational case, the orbital distance de- 
creases gradually in the initial stages, and then when the 
orbital instability is triggered, it quickly decreases lead- 
ing to merger. However, the behavior of the merger is 
different from that in the corotational cases. For the 
irrotational binary, the initial distribution of angular ve- 
locity around the center of mass is a decreasing function 



of the distance from the center (and the absolute value 
of the velocity is almost independent of position; cf. 
Figs. 9-11 at t = 0). Hence, the centrifugal force in the 
outer region of the merged objects is not as strong as 
that in the corotational cases. Consequently, spiral arms 
are not formed in a significant way. On the other hand, 
the magnitude of the centrifugal force in the inner region 
is stronger than that in the corotational cases. As a re- 
sult, two oscillating cores are formed in the inner region, 
and this structure is maintained for a short while. These 
features have been found also in Newtonian simulations 

In the case of model (II), the two cores bounce after 
their first collision, and then they merge to form an oscil- 
lating new neutron star. In Figs. 12 and 13, we show the 
density contour lines for in the x — z plane and the an- 
gular velocity defined hy fl = {xv^ — yv^ ) / {x^ -f ) along 
the X and y-axes in the equatorial plane at t = 1.81Porb- 
We find that the new neutron star has a toroidal struc- 
ture which is sustained by differential rotation [|T] . Note 
that 99.5% of the total rest mass is inside r ~ 6A/go and 
appears to constitute the merged object at t = 1.81Poib 
in this case. Thus, the rest mass of the new neutron star 
is ~ 45% larger than Af^^^ax As in the corotating 
case, K' increases only by a small factor (at most 10%) 
in the high density region, so that the role of the ther- 
mal pressure increase for supporting the large mass is not 
significant. 

In this model (II), the initial value of J/Af^ is less than 
unity and the final value should be even smaller as argued 
in Sec. IV. A. Nevertheless, the merged object does not 
collapse to form a black hole. Axisymmetric simulations 
of stellar core collapse |Q have indicated that a black 
hole is formed for J/Mg < 1 in most cases if the mass is 
large enough, and that the angular momentum parameter 
is a good indicator for predicting the final product. The 
present simulation suggests that this is not always the 
case for the merger of binary neutron stars. 

We note again that the new rotating neutron star 
was non-axisymmetric when we stopped the simulation. 
Therefore it will evolve secularly and may become unsta- 
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ble against gravitational collapse to a black hole after a 
substantial amount of the angular momentum has been 
carried away by gravitational radiation ||37|] . 

For model (13), the inner part contracts due to self- 
gravity without bouncing because the pressure and cen- 
trifugal force are not strong enough to balance the self- 
gravity. Consequently, it quickly collapses to form a black 
hole. For model (12), on the other hand, the self-gravity 
is weaker than that for model (13), so that the two cores 
bounce at the first collision for t ^ 1.2Poib- Then, they 
approach again redistributing the angular momentum as 
well as dissipating it by gravitational radiation, and fi- 
nally the merged object forms a black hole. To demon- 
strate this feature, we show a at r = as a function of 
t/Povh in Fig. 14. For model (13), a{r = 0) monoton- 
ically approaches zero, but for model (12), it increases 
again after it reaches a first minimum at t ^ 1.2Porb- 
These numerical results indicate that the merging pro- 
cess towards the final state depends considerably on the 
initial compactness of the neutron stars. 

In Fig. 15, we show the time evolution of the fraction 
of the rest mass inside a coordinate radius r, Af,(r)/Af,, 
for models (12) and (13). We again choose r — 1.5, 3 and 
4.5Mgo as coordinate radii. It is found that more than 
99% of the total rest mass was inside r = 4.5Afgo for 
both models when we stopped the simulations. Thus, in 
contrast with the corotational cases, only a tiny fraction 
of the total rest mass (< 1%) can form a disk around 
the black hole at r > 4.5Mgo- For model (13), > 99% 
of the total rest mass is inside r = 1.5Afgo which almost 
coincides with the location of the apparent horizon at 
the final snapshot of Fig. 11. Hence, we can conclude 
that the disk mass is very small (< O.OlAf*) for model 
(13). For model (12), we could not determine the location 
of the apparent horizon before the simulation crashed, 
and so we cannot make any strong conclusion. However, 
Fig. 15 shows that the mass fraction outside r — 3AIgo 
is < O.OlAf* and that inside r = l.dMgo is quickly in- 
creasing at i ~ l.SPoib- Hence, the final disk mass again 
appears to be very small as in the case (13). 



C. gravitational waves 

To extract gravitational waveforms, we define non- 
dimensional variables 

h+ = r(7,, - jyy)/{2Mgo), (4.4) 
hx = rjxy/Mgn, (4.5) 

along the z-axis. Since we adopt the AMD gauge con- 
dition and have prepared initial conditions for which 
^^•'dijjk = 0, "fij is approximately transverse and trace- 
less in the wave zone [E7|. As a result, /i+ and h X arc ex- 
pected to be appropriate measures of gravitational waves. 

In Figs. 16, we show waveforms for corotational models 
(CI) (the solid lines) and (C2) (the dashed lines), and in 
Figs. 17, for irrotational models (II) (the solid lines) and 
(12) (the dashed lines) as a function of retarded time 
{t — 2;obs)/-Porb where Zobs denotes the point along the 
z-axis at which the waveforms are extracted. 

For both corotational and irrotational cases, the am- 
plitude gradually rises with decreasing orbital separation, 
but after the amplitude reaches the maximum, the wave- 
forms for the two cases have different characters. In the 
corotational cases, the amplitude soon becomes small af- 
ter the maximum, while in the irrotational cases, it does 
not become small very quickly, but has a couple of fairly 
large peaks. The reason is that the double core structure 
which enhances the amplitude is preserved for a short 
while after the merger in the irrotational cases. Such 
a feature has been found also in Newtonian simulations 
||l^,0, indicating that the Newtonian simulations are 
helpful for investigation of the qualitative outcome of 
gravitational waveforms. In particular, the waveforms 
for models (CI) and (II) in which new neutron stars 
are formed are qualitatively similar to the corresponding 
Newtonian models [ p^|Jl^ , although quantitative features 
such as amplitude and wavelength are different. There- 
fore, the Newtonian simulation is useful as a guideline for 
fully GR simulations particularly when the final product 
is a neutron star. 

The maximum amplitude for /i+ and hx is typically 
0.1 as shown in Figs. 16 and 17. This implies that the 
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typical maximum amplitude of gravitational waves from 
a source at the distance r is 

As we mentioned above, the outer boundaries of the 
computational domain with 233x233x117 grid points are 
located at ^ O.SAgw on each axis. This imphes that the 
waveforms extracted are not accurate asymptotic wave- 
forms. For example, a slight unrealistic modulation (the 
wave amplitude deviates gradually with time in the pos- 
itive direction) is found in the waveform for in every 
case which seems due to numerical error. 

To estimate the magnitude of the error, we performed 
one large simulation for model (12) with grid size 293 x 
293 X 147, fixing the grid spacing but widening the com- 
putational region. Even in this case, the outer boundaries 
are located at ~ 0.35Agw on each axis. In Fig. 18, we 
show the waveforms for 293 x 293 x 147 (the solid lines) , 
233 X 233 X 117 (the dashed lines), and 193 x 193 x 97 
(the dotted lines). For the early phase of merging, the 
magnitude of the modulation is smaller and smaller with 
increasing number of grid points, which implies that this 
effect is spurious due to the restricted computational re- 
gion. For the very late phase of merging, on the other 
hand, the magnitude of the modulation does not change 
even with widening the computational domain. This sug- 
gests that the resolution of the central regions of the 
merged object is not sufficient in that phase to compute 
accurate waveforms. We also find that the wave ampli- 
tude increases slightly with widening the computational 
region. This indicates that the amplitudes shown in Figs. 
16 and 17 might underestimate the asymptotic one by 
several tens of percent especially in the early phase of 
the merger. All of these facts indicate that we need a 
larger scale computation to improve the accuracy of the 
gravitational waveforms. 

Even in the case of black hole formation, the shapes of 
the waveforms are similar to those in the neutron star for- 
mation case before the gravitational collapse to a black 
hole has occurred (compare the waveforms for models 
(II) and (12)). The difference in the waveforms will ap- 



pear after the gravitational collapse. However, since we 
could not continue simulations for a long time at this 
stage, we cannot describe the features of the waveform 
in detail. In the following, we speculate on the expected 
outcome and discuss the significance of the waveforms 
from the observational point of view. 

According to the standard scenario, the quasi normal 
modes of the black hole are excited in the final phase of 
black hole formation, and the amplitudes of these modes 
subsequently damps. As we showed in the previous two 
subsections, the formation timescale of the black hole 
is different depending on the compactness of the neu- 
tron stars before the merger. This implies that the time 
duration from the moment when the amplitude of the 
gravitational waves becomes a maximum to the moment 
when the amplitude of the waves from the merged object 
damps depends on the initial compactness parameter of 
the neutron stars (see Fig. 19, which shows a schematic 
picture for the expected gravitational waveforms) : In the 
case of neutron star formation (c/. Fig. 19 (a)), the 
damping time for quasi-periodic gravitational waves of 
small amplitude emitted from non-axisymmetric defor- 
mation of the new neutron star is the timescale of gravi- 
tational radiation reaction which is much longer than the 
dynamical timescale. In the case of black hole formation, 
we have a number of possibilities: If the compactness pa- 
rameter of the neutron stars before the merger is not very 
large (c/. Fig. 19 (b)), the timescale for the formation 
process is fairly long and the quasi-periodic oscillations 
due to non-axisymmetric deformation of the merged ob- 
ject will be seen for a short while after the merger. If 
the neutron stars are sufficiently compact (c/. Fig. 19 
(c)), the black hole is formed quickly and the amphtude 
of gravitational waves will also damp quickly. Therefore, 
the time duration from the gravitational wave burst to 
its damping (note that we do not need here the detail 
of the waveforms) will constrain the initial compactness 
of the neutron stars, and, consequently, the equation of 
state for high density neutron matter 
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V. DISCUSSION 

As we found in Sec. IV, the final products of the 
merger depend sensitively on the initial compactness of 
the neutron stars. In the corotational case, (1) the final 
product is a massive neutron star when the ratio of the 
rest mass of each star to Ml^max (C'mass) is 0.8 ; (2) 
the final product is a black hole when Cmass is ^ 0.9. 
If it is at most 0.9, the formation timescale is longer 
than the dynamical timescale (or the oscillation period 
of the merged object, Pose)- On the other hand, if Cmass 
is ~ 1, the formation timescale of the black hole is as 
short as the dynamical timescale Pose)- In the irro- 
tational case, (3) the final product is a massive neutron 
star when Cmass is ^ 0.7 ; (4) the final product is a black 
hole when Cmass is 0.8. If it is at most ~ 0.8, the for- 
mation timescale is longer than Pose- On the other hand, 
if Cmass is larger than ^ 0.9, the formation timescale is 
< P 

^ OSC ■ 

Let us consider the case where two irrotational neutron 
stars of rest mass 1.6Mq (i.e., the gravitational mass 
is 1.4Mq) merge. The numerical results in this pa- 
per indicate that (a) if Af^^max is less than ^ 1.8Mq, 
the merged object forms a black hole on the dynamical 
timescale Pose ; (b) if Af^^max is ^ 2A^q, the final 
product is also a black hole, but the formation timescale 
is longer than Pose; (c) if A/^'^max is larger than ^ 2.2AfQ, 
the final product will be a massive neutron star. This fact 
provides us the following interesting possibility. Suppose 
that we will be able to find the mass of each neutron star 
during the inspiraling phase by means of the matched 
filtering method |Q with the aid of the post-Newtonian 
template Then, if we observe the merger process 

to the final products, in particular the timescale for for- 
mation of a black hole, we can constrain the maximum 
allowed neutron star mass, and consequently the nuclear 
equation of state. 

Unfortunately, the frequency of gravitational waves af- 
ter the merger will be so high (typically Po^^;^ ^ ^-^mh ~ 
2 — SkHz) that laser interferometers such as LIGO ||] 
will not be able to detect them. To observe such high 



frequency gravitational waves, specially designed nar- 
row band interferometers or resonant-mass detectors are 
needed ||]. We should keep in mind that such future 
gravitational wave detectors would have the possibility 
to provide us with important information about the neu- 
tron star equation of state. 

Another important outcome of the present simulations 
concerns the mass of the disk around a black hole formed 
after the merger. A disk of mass ^ 0.05 — O.IA/* may be 
formed around a black hole after the merger of corota- 
tional binaries. However, for the merger of irrotational 
binaries, the mass of the disk appears to be very small 
< O.OlAf*. An irrotational velocity field is considered 
to be a good approximation for realistic binary neutron 
stars before merger Therefore, a massive disk of 

mass > O.IAz/q may not be formed around a black hole 
after the merger of binary neutron stars of nearly equal 
mass. This is not very promising for some scenarios for 
GRBs, in which a black hole-toroid system formed after 
the merger of nearly equal mass binary neutron stars is 
considered to be its central engine. 

We have performed simulations using a modified form 
of the ADM formalism for the Einstein field equation 
with the AMD gauge and approximate maximal slicing 
conditions Needless to say, simulations by other 

groups using different formulations, gauge conditions and 
numerical implementations ||l9|j20|| are necessary to re- 
confirm the present results. 

In this paper, we have performed simulations only for 
the case F = 2. As Newtonian simulations have indicated 
p^ , p^ , the merging process and final products may also 
depend sensitively on the stiffness of neutron star mat- 
ter. In a forthcoming paper, we will perform simulations 
changing T to investigate the dependence on the stiffness 
of the equation of state and to clarify whether the present 
conclusions are modified or not. 
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Table I. A list of several quantities for initial conditions of binary neutron stars. The maximum density, total 
rest mass M*, gravitational mass Mgo, J/M^g, compactness Cj = (M(,ofi)^^'^(~ Mgo/a where a is orbital separation), 
ratio of the emission timescale of gravitational waves to the orbital period Rt = 5{Mgo^)~^^"^ / 128tt , the ratio of the 
rest mass of each star to the maximum allowed mass for a spherical star Cmass = Af*/2M*'^max, the ratio of Mgo to 
the grid spacing in the simulation {Mgo/ Ax), the type of velocity field and final products are shown. Here, Ml^'max 
denotes the maximum allowed mass for a spherical star (~ 1.435). Here, we quote values for K = 200/7r. The mass 
and density can be scaled by the rules M*(ii'7r/200)^/2, Mgo{K'!r/200y/'^, and Pn,^^{KTr/200)-'^ for arbitrary K, while 
other non-dimensional quantities are invariant. 
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* We initially reduced the angular momentum from the corresponding quasi-equilibrium states by ~ 2.5%. 




FIG. 1. Rest mass M* as a function of maximum density pmax for each star in the binary for K = 200/7r. The binaries which 
are used in the simulations are marked with (CI), (C2), (C3), (II), (12), and (13). The solid line denotes the relation for the 
spherical stars. We note that the mass and the density can be scaled by the rules M*(i^7r/200)^/^ and pmax(A"7r/200)"^ for 
arbitrary K. Scales for the top and right axes are shown for K = 123.7 in which the maximum rest mass for spherical stars is 
2M0. 
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FIG. 2. Snapshots of the density contour lines for p, and the velocity field {v^, v^) in the equatorial plane for model (CI). 
The contour lines axe drawn for p*/p^ max = lO"'''^-', where p, max denotes the maximum value of p, at t = (here it is 
0.00441), for J = 0, 1, 2, • • • , 10. Vectors indicate the local velocity field and the scale is as shown in the top left-hand frame. P 
denotes the initial orbital period Porb. The length scale is shown in units of GMgo/c? . 
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FIG. 3. The same as Fig. 2, but for model (C2). The contour lines are drawn for p«/p* 




10 



-0.3j 



where 



p* max = 0.00757, for j = 0, 1, 2, • • • , 10. The dashed line in the last figure denotes the circle with r = A.5Mgo within 
which ~ 95% of the total rest mass is included. 
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FIG. 4. The same as Fig. 2, but for model (C3). The contour lines are drawn for p«/p* max = 10~°'^^ , where p* max = 0.0171, 
for J = 0, 1, 2, ■ • • , 10. The dashed line in the last snapshot denotes the circle with r = 3Mgo within which ~ 95% of the total 
rest mass is included. The thick solid line for r ~ Mgo in the last snapshot denotes the location of the apparent horizon. Note 
that there are ~ 7 grid points along the radius of the apparent horizon. 
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FIG. 5. The density contour lines for p* in the y = plane at t = 2.07Porb for model (CI). The contour lines are drawn in 
the same way as in Fig. 2. The length scale is shown in units of GMgo/(? . 
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FIG. 6. The angular velocity Q. along the x-axis (solid line) and y-axis (dotted line) at f = 2.07Porb for model (CI). The 
length scale and f2 are shown in units of GMgo/(? and <? /GMgo, respectively. 
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FIG. 7. Q at r = as a function of t/Porb for models (CI), (C2) and (C3). 




FIG. 8. Fraction of the rest mass inside a coordinate radius r as a function of t/Porh for models (C2) and (03) in which a 
black hole is formed after the merger. 
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FIG. 10. The same as Fig. 2, but for model (12). The contour lines are drawn for p«/p* max = lO"" '^-' , where 
p* max = 0.00642, for j = 0, 1, 2, • • • , 10. The dashed line in the last figure denotes the circle with r = A.hMgo within 
which more than 99% of the total rest mass is included. 
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FIG. 11. The same as Fig. 2, but for model (13). The contour lines are drawn for p«/p* max ~ lO^" "'-' , where 
p* max = 0.0136, for j = 0, 1, 2, • • • , 10. The dashed line in the last snapshot denotes the circle with r = 3Mgo within 
which more than 99% of the total rest mass is included. The thick solid line for r ~ Mgo in the last snapshot denotes the 
location of the apparent horizon. Note that there are ~ 7 grid points along the radius of the apparent horizon. 
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FIG. 12. The density contour lines for p« in the y = plane at t = l.SlPorb for model (II). The contour lines are drawn in 
the same way as in Fig. 9. The length scale is shown in units of GMgo/(? . 




FIG. 13. The angular velocity Q. along the x-axis (solid line) and t/-axis (dotted line) at f = l.SlPorb for model (II). The 
length scale and Q, are shown in units of GMgo/(? and <? /GMgo, respectively. 
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FIG. 14. a at r = as a function of t/Pah for models (II), (12) and (13). 




FIG. 15. Fraction of the rest mass inside a coordinate radius r as a function of t/Poib for models (12) and (13) in which a 
black hole is formed after the merger 
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FIG. 17. h+ and hx as functions of retarded time for the irrotational models (II) (solid lines) and (12) (dashed lines). 
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FIG. 18. /i+ and as functions of retarded time for irrotational models (12) with 293 x 293 x 147 grid size (solid lines), 
233 X 233 X 117 grid size (dashed lines), and 193 x 193 x 97 grid size (dotted lines). In each case, the outer boundaries (and 
the points where the waveforms are extracted) are located at z ~ 0.35, 0.28 and 0.23Agw, respectively. 




FIG. 19. Schematic pictures for expected gravitational waveforms during and after the merger for (a) the neutron star 
formation case; (b) the black hole formation case in which the compactness of the neutron stars before the merger is not very 
large and the formation timescale is fairly long; (c) the black hole formation case in which the compactness of the neutron stars 
before the merger is large enough that the formation timescale is short. 
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